/* -----------------------------------------------------------------------
   See COPYRIGHT.TXT and LICENSE.TXT for copyright and license information
   ----------------------------------------------------------------------- */
#include "plm_config.h"
#include <math.h>
#include <string.h>
#include <stdlib.h>
#include "ramp_filter.h"
#include "fftw3.h"
#define MARGIN5 5

#ifndef PI
static const double PI = 3.14159265;
#endif

#ifndef DEGTORAD
static const double DEGTORAD = 3.14159265 / 180.0;
#endif

#ifndef MARGIN
static const unsigned int MARGIN = 5;
#else
#error "MARGIN IS DEFINED"
#endif

/* GCS: Dec 12, 2009 
   Change from unsigned short to float,
   Change to destructive update
*/
//! In-place ramp filter for greyscale images
//! \param data The pixel data of the image
//! \param width The width of the image
//! \param height The height of the image
//! \template_param T The type of pixel in the image
void
ramp_filter(
    float *data, 
    unsigned int width, 
    unsigned int height
)
{
    unsigned int i, r, c;
    unsigned int N;

    fftw_complex *in;
    fftw_complex *fft;
    fftw_complex *ifft;
    fftw_plan fftp;
    fftw_plan ifftp;
    double *ramp;
    ramp = (double*) malloc (width * sizeof(double));
    if (!ramp) {
        print_and_exit ("Error allocating memory for ramp\n");
    }
    N = width * height;
    in = (fftw_complex*) fftw_malloc (sizeof(fftw_complex) * N);
    fft = (fftw_complex*) fftw_malloc (sizeof(fftw_complex) * N);
    ifft = (fftw_complex*) fftw_malloc (sizeof(fftw_complex) * N);
    if (!in || !fft || !ifft) {
        print_and_exit ("Error allocating memory for fft\n");
    }

    for (r = 0; r < MARGIN; ++r)
        memcpy (data + r * width, data + MARGIN * width, 
	    width * sizeof(float));

    for (r = height - MARGIN; r < height; ++r)
        memcpy (data + r * width, data + (height - MARGIN - 1) * width, 
	    width * sizeof(float));

    for (r = 0; r < height; ++r) {
        for (c = 0; c < MARGIN; ++c)
            data[r * width + c] = data[r * width + MARGIN];
        for (c = width - MARGIN; c < width; ++c)
            data[r * width + c] = data[r * width + width - MARGIN - 1];
    }

    // Fill in
    for (i = 0; i < N; ++i) {
        in[i][0] = (double)(data[i]);
        //in[i][0] /= 65535;
        //in[i][0] = (in[i][0] == 0 ? 1 : in[i][0]);
        //in[i][0] = -log (in[i][0]);
        in[i][1] = 0.0;
    }

    for (i = 0; i < width / 2; ++i)
        ramp[i] = i;

    for (i = width / 2; i < (int) width; ++i)
        ramp[i] = width - i;

    for (i = 0; i < width; ++i)
        ramp[i] *= (cos (i * DEGTORAD * 360 / width) + 1) / 2;

    for (r = 0; r < height; ++r)
    {
	fftp = fftw_plan_dft_1d (width, in + r * width, fft + r * width, 
	    FFTW_FORWARD, FFTW_ESTIMATE);
	if (!fftp) {
	    print_and_exit ("Error creating fft plan\n");
	}
	ifftp = fftw_plan_dft_1d (width, fft + r * width, ifft + r * width, 
	    FFTW_BACKWARD, FFTW_ESTIMATE);
	if (!ifftp) {
	    print_and_exit ("Error creating ifft plan\n");
	}

        fftw_execute (fftp);

        // Apply ramp
        for (c = 0; c < width; ++c) {
            fft[r * width + c][0] *= ramp[c];
            fft[r * width + c][1] *= ramp[c];
        }

        fftw_execute (ifftp);

	fftw_destroy_plan (fftp);
	fftw_destroy_plan (ifftp);
    }

    for (i = 0; i < N; ++i)
        ifft[i][0] /= width;

    for (i = 0; i < N; ++i)
        data[i] = (float)(ifft[i][0]);

    fftw_free (in);
    fftw_free (fft);
    fftw_free (ifft);
    free (ramp);
}



	void RampFilter_filter_extend_DC(float * data,
		float* out,
		unsigned int prewidth,
		unsigned int height,
		int nproj)
	{	unsigned int i,j,r,c;
	static unsigned int N;
	fftw_complex *in;
	fftw_complex *fft; 
	fftw_complex *ifft;
	fftw_plan fftp;
	fftw_plan ifftp;
	double* ramp;
	unsigned int width=512;
	unsigned head=ceil(0.67*(width-prewidth));
	unsigned tail=width-head;

	ramp=(double*)malloc(width*sizeof(double));
	if (ramp==NULL){
		printf("Malloca error");
		exit(1);
	}

	N= width * height;
	in    = (fftw_complex*)fftw_malloc(sizeof(fftw_complex) * N);
	fft   = (fftw_complex*)fftw_malloc(sizeof(fftw_complex) * N);
	ifft  = (fftw_complex*)fftw_malloc(sizeof(fftw_complex) * N);

	//for (r=0; r<MARGIN5; ++r)
	//	memcpy(data+r*width,data+MARGIN5*width,width*sizeof(unsigned short ));

	//for (r=height-MARGIN5; r<height; ++r)
	//	memcpy(data+r*width,data+(height-MARGIN5-1)*width,width*sizeof(unsigned short ));

	//for (r=0; r<height; ++r){
	//	for (c=0 ; c<MARGIN5; ++c)
	//		data[r*width+c]=data[r*width+MARGIN5];
	//	for (c=width-MARGIN5 ; c<width; ++c)
	//		data[r*width+c]=data[r*width+width-MARGIN5-1];
	//}  

	// Fill in
	//for (i = 0; i < N; ++i)
	//{
	//	in[i][0] = (double)(data[i]);
	//	in[i][1] = 0.0;
	//}

	for (j = 0; j < height; j++)
	{
		for (i=0; i<head; i++){
			in[width*j+i][0] = (double)(data[prewidth*j+0]);
			in[i][1] = 0.0;
		}
		for (i=head; i<head+prewidth; i++){
			in[width*j+i][0] = (double)(data[prewidth*j+i-head]);
			in[i][1] = 0.0;
		}
		for (i=head+prewidth; i<width; i++){
			in[width*j+i][0] = (double)(data[prewidth*j+prewidth-1]);
			in[i][1] = 0.0;
		}
	}



	for (i = 0; i < width / 2; ++i)
		ramp[i] = i;
	//ramp[0]=100.0f/((float)nproj);

	for (i = width / 2; i < (int)width; ++i)
		ramp[i] = width - i;

	for (i = 0; i < width; ++i)
		ramp[i] *= (cos(i * DEGTORAD * 360/width)+1)/2;

	for (r = 0; r < height; ++r)
	{
		fftp  = fftw_plan_dft_1d(width, in + r*width, fft + r*width, FFTW_FORWARD, FFTW_ESTIMATE);
		ifftp = fftw_plan_dft_1d(width, fft + r*width, ifft + r*width, FFTW_BACKWARD, FFTW_ESTIMATE);

		fftw_execute(fftp);

		// Apply ramp
#if 1
		for (c = 0; c < width; ++c)
		{
			fft[r*width + c][0] *= ramp[c];
			fft[r*width + c][1] *= ramp[c];
		}
#endif

		fftw_execute(ifftp);

		fftw_destroy_plan(fftp);
		fftw_destroy_plan(ifftp);
	}


	for (i = 0; i < N; ++i)
		ifft[i][0] /= width;

	//for (i = 0; i < N; ++i)
	//	out[i] = (float)(ifft[i][0]);

	for (j = 0; j < height; j++)
		for (i=head; i<head+prewidth; i++)
			out[prewidth*j+i-head] = (double)(ifft[width*j+i][0]);



	fftw_free(in);
	fftw_free(fft);
	fftw_free(ifft);
	free(ramp);
	}



